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We perform fully general relativistic simulations of rotating stellar core collapse in three spatial 
dimension. The hydrodynamic equations are solved using a high-resolution shock-capturing scheme. 
A parametric equation of state is adopted to model collapsing stellar cores and neutron stars following 
Dimmelmeier et al. The early stage of the collapse is followed by an axisymmetric code. When the 
stellar core becomes compact enough, we start a three-dimensional simulation adding a bar-mode 
nonaxisymmetric density perturbation. The axisymmetric simulations are performed for a wide 
variety of initial conditions changing the rotational velocity profile, parameters of the equations of 
state, and the total mass. It is clarified that the maximum density, the maximum value of the 
compactness, and the maximum value of the ratio of the kinetic energy T to the gravitational 
potential energy W {(3 = T/W) achieved during the stellar collapse and bounce depend sensitively 
on the velocity profile and the total mass of the initial core, and equations of state. It is also 
found that for all the models with high degree of difi'erential rotation, a funnel structure is formed 
around the rotational ax;is after the formation of neutron stars. For selected models in which the 
maximum value of (3 is larger than ~ 0.27, three-dimensional numerical simulations are performed. 
It is found that the bar-mode dynamical instability sets in for the case that the following conditions 
are satisfied: (i) the progenitor of the stellar core collapse should be rapidly rotating with the initial 
value of 0.01 ^ /3 ^ 0.02, (ii) the degree of differential rotation for the velocity profile of the initial 
condition should be sufficiently high, and (iii) a depletion factor of pressure in an early stage of 
collapse should be large enough to induce a significant contraction to form a compact stellar core for 
which an efficient spin-up can be achieved surmounting the strong centrifugal force. As a result of 
the onset of the bar-mode dynamical instabilities, the amplitude of gravitational waves can be by a 
factor of ~ 10 larger than that in the axisymmetric collapse. It is found that a dynamical instability 
with the m — 1 mode is also induced for the dynamically unstable cases against the bar-mode, but 
the perturbation does not grow significantly and, hence, it does not contribute to an outstanding 
amplification of gravitational waves. No evidence for fragmentation of the protoneutron stars is 
found in the first a few 10 msec after the bounce. 
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I. INTRODUCTION 

One of the most important issues of hydrodynamic simulations in general relativity is to clarify stellar core collapse 
to a neutron star or a black hole. The formation of neutron stars and black holes is among the most promising 
sources of gravitational waves. This fact has stimulated numerical simulations for the stellar core collapse [1-12]. 
Hovirever, most of these works have been done in the Newtonian framework and in the assumption of axial symmetry. 
As demonstrated in [10,12], general relativistic effects modify the dynamics of the collapse and the gravitational 
waveforms significantly in the formation of neutron stars. Thus, the simulation should be performed in the framework 
of general relativity. The assumption of axial symmetry is appropriate for the case that the rotating stellar core is not 
rapidly rotating. However, for the sufhciently rapidly rotating cases, nonaxisymmetric instabilities may grow during 
the collapse and the bounce [7]. As a result, the amplitude of gravitational waves may be increased significantly. 

To date, there has been no general relativistic work for the stellar core collapse in three spatial dimensions. Three- 
dimensional simulations of the stellar core collapse have been performed only in the framework of Newtonian gravity 
[4,7]. Hydrodynamic simulations for gravitational collapse or for the onset of nonaxisymmetric instabilities of rotating 
neutron stars in full general relativity have been performed so far [13-17], but no simulation has been done for the 
rotating stellar core collapse to a neutron star or a black hole. In this paper, we present the first numerical results of 
three-dimensional simulations for rapidly rotating stellar core collapse in full general relativity. 

Three-dimensional simulation is motivated by two major purposes. One is to clarify the criterion for the onset 
of nonaxisymmetric dynamical instabilities during the collapse, and the outcome after the onset of the instabilities. 
So far, a number of numerical simulations have illustrated that rapidly rotating stars in isolation and in equilibrium 
are often subject to nonaxisymmetric dynamical instabilities not only in Newtonian theory [18-28], but also in post- 
Newtonian approximation [29], and in general relativity [15]. These simulations have shown that the dynamical 
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bar-mode instabilities set in (i) when the ratio of the kinetic energy T to the gravitational potential energy W 
(hereafter [3 = T/W) is larger than ~ 0.27 or (ii) when the rotating star is highly differentially rotating, even for 
(3 <C 0.27 [28]. As a result of the onset of the nonaxisymmetric instabilities, a bar and spiral arms are formed 
which can redistribute angular momentum profile and change the density profile of the star. Also, a burst-type and 
subsequent quasiperiodic gravitational waves with a high amplitude can be emitted in the case of rapidly rotating 
neutron stars [22,28,15]. However, the numerical simulations have been performed mostly for isolated rotating stars in 
equilibrium. To our knowledge, [7] is only one published paper in which the nonaxisymmetric dynamical instabilities 
during stellar core collapse have been investigated. In [7], the authors performed Newtonian simulations for a few 
models and indicated that the dynamical instability sets in only for the case that the value of /3 exceeds much beyond 
0.27. Such condition is satisfied only when the progenitor of the collapse is rapidly and highly differentially rotating 
and the depletion of the internal energy in an early stage of the collapse is large enough to produce a very compact 
core for which a significant spin-up can be achieved surmounting the strong centrifugal force [6]. 

Although the previous Newtonian work [7] indicated a criterion for the onset of nonaxisymmetric dynamical in- 
stabilities, many unclear points still remain unsolved as follows. First, in the Newtonian analysis [7], the authors 
adopted a parametric equation of state, and performed simulations changing its own parameters. They found that for 
the onset of the nonaxisymmetric instability during collapse, a soft equation of state with Fi = 1.28 and F2 = 2.5 is 
necessary (see Sec. II B for the definition of Fi and r2). Unfortunately, in the equations of state that they adopted, 
the maximum gravitational mass for cold spherical neutron stars in general relativity becomes « I.SMq (see Table I), 
which is too small to be adopted as a plausible equation of state in general relativistic simulations since the maximum 
mass of neutron stars for a given equation of state should be larger than « I.MMq which is the precisely determined 
mass of a neutron star in PSRB1913-I-16 [30]. A study with more plausible equations of state is required. 

Second, the authors in [7] focus little on the instabilities associated with m = 2 bar mode, although it is the fastest 
growing mode of the nonaxisymmetric dynamical instabilities for equilibrium stars in most cases. (Here, m denotes 
the azimuthal quantum number.) Thus, the criterion for the onset of the bar-mode instabilities are not still clear. 
Also, they paid little attention to the bar-mode dynamical instabilities for highly differentially rotating cases such as 
those recently reported in [28]. This instability can set in even for a small value of /3 < 0.27. This implies that for 
highly differentially rotating initial conditions, attention should be also paid for small values of f3. 

Third, in general relativity, the collapsed core can reach a more compact state than that simulated in the Newtonian 
theory due to the fact that self-gravity becomes stronger [10]. As a result, more efficient spin-up will be achieved. 
Therefore, the probability for the onset of nonaxisymmetric dynamical instabilities would be underestimated in the 
Newtonian simulation. This suggests that general relativistic analysis may be crucial for the study of nonaxisymmetric 
dynamical instabilities. 

Finally, in [7], the mass of the stellar core adopted is set to be in a narrow range between 1.5 and 1.7M0. According 
to the theory of stellar evolutions, in a very massive star of low metallicity with the initial mass 5OM0 ^ M ^ lOOAf©, 
the produced iron core may become 2-3Mq [31-33]. This indicates that the mass of the core in nature may be in a 
wide range between ^ IMq and ~ 3Mq. With the increase of the mass, the self-gravity becomes stronger, and hence, 
the collapsed stellar core can reach a more compact state for which a spin-up may be enhanced effectively. Thus, the 
larger core mass may increase the probability for the onset of nonaxisymmetric dynamical instabilities. 

Motivated by the questions mentioned above, we perform general relativistic simulations choosing rapidly and 
highly differentially rotating massive stars with plausible equations of state and with a wide mass range. Following 
[7], we adopt a parametric equation of state. However, we choose sets of the parameters in which the maximum 
Arnowitt-Descr-Misner (ADM) mass of cold spherical neutron star becomes w 1.6Mq. With this setting, we choose 
the mass of the stellar core in the range between ~ 1.5 and ~ 3Mq. Furthermore, we pay particular attention to the 
bar-mode instabilities which are likely to be the fastest growing mode. 

Another major role of three-dimensional simulations for the stellar core collapse is to determine the amplitude 
and the characteristic frequency of gravitational waves in the case that the nonaxisymmetric dynamical instabilities 
set in. Axisymmetric numerical simulations have clarified that the amplitude of gravitational waves emitted in the 
stellar core collapse is at most several xlO~^'^ at a distance of 10 Mpc (e.g., [10]), and the frequency is between 100 
Hz and 1 kHz. Although the frequency is in the most sensitive band of the laser interferometric gravitational wave 
detectors [34] . the value of the amplitude is too small to be detected if the stellar core collapse occurs outside our 
local group of galaxies. In the three-dimensional process, on the other hand, the amplitude is often by a factor of 
~ 10 larger than that in the axisymmetric phenomena because of the increase of the degree of asymmetry. Hence, if 
the nonaxisymmetric instabilities set in, the stellar core collapse may become a much stronger emitter of gravitational 
waves than that considered so far. 

This paper is organized as follows. In Sec. II, we briefly review our formulation of general relativistic simulation, 
equations of state adopted in this paper, and methods for extraction of gravitational waves. In Sec. Ill, initial 
conditions and computational setting are described. In Sec. IV, numerical results of axisymmetric simulations are 
presented paying attention to the value of (3 and to the profiles of the density and the angular velocity of the outcomes. 
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In Sec. V, numerical results of three-dimensional simulations are presented, clarifying the criterion for the onset of bar- 
mode dynamical instabilities. Gravitational waveforms emitted in the growth of the bar- mode dynamical instabilities 
are also presented. Sec. VI is devoted to a summary. Throughout this paper, we adopt the geometrical units in which 
G = c = 1 where G and c are the gravitational constant and speed of light, respectively. 



II. FORMULATION 



A. Summciry of basic equations and implementations 

We perform hydrodynamic simulations in full general relativity using the same formulation as in [35,36], to which 
the reader may refer for details and basic equations. The fundamental variables for the hydrodynamics are p : rest 
mass density, e : specific internal energy, P : pressure, : four velocity, and 

where subscripts k, - ■ ■ denote x, y and z, and /z the spacetime components. As the fundamental variables to be 
evolved in the numerical simulations, we define a weighted density p*, a weighted four- velocity Wj, and a specific 
energy density e as 

= ptwe^^, (2) 
Ui = hui, (3) 
P 

e = hw , (4) 

pw 

where e"^ denotes the conformal factor, w = au* , and h= 1 + e + P/ p. Here, e is computed from T^^n'^n'^ /{pw) where 
Tfj_i, and n'' denote the energy-momentum tensor and a timelike unit normal vector. General relativistic hydrodynamic 

equations are solved using a high-resolution shock-capturing scheme [37,35]. In axisymmetric and three-dimensional 
simulations, the cylindrical and Cartesian coordinates are used, respectively. The details of our hydrodynamic code 
arc described in [35]. 

For the following, we define the total baryon rest-mass, internal energy, and rotational kinetic energy of the system 

as 



M, 
U 
T 



* = J d^xp^ , (5) 
= j (fxp*e, (6) 
= ^ j (fxp^u^v'^, (7) 



where is the conserved quantity. The definitions of U and T agree with those for axisymmetric rotating stars in 
equilibrium [39] . In the axisymmetric case, the angular momentum J is a conserved quantity, and defined by 



J 



= J (fxp^u^. (8) 



Note that in the nonaxisymmetric case, the equation of J has a different form (e.g., [13]). 

The fundamental variables for the geometry are a: lapse function, /S'': shift vector, jij: three-metric, 7 = e^^"^ = 
det(7y): trace of the three-metric, jij = e^'^'^jij: conformal three-metric, and Kij : extrinsic curvature. We evolve 
7ij, (j), Aij = e~^'^{Kij — ^ijKjJ'), and the trace of the extrinsic curvature K/} together with the three auxiliary 
variables Fi = S^'^dj'jik with an unconstrained free evolution code as in [40,13,41,35]. The Einstein equations are 
solved in the Cartesian coordinates. In the axisymmetric case, the Cartoon method is used [42,43]. In both cases, the 
equatorial reflection symmetry is assumed. The outer boundary conditions we adopt are the same as in the previous 
papers (e.g., [13,41,35]). 

As the slicing condition, we impose an "approximate" maximal slice condition (KfJ' w 0) which is the same as that 
adopted in previous papers (e.g., [13,15,41]). As the spatial gauge condition, we adopt a hyperbolic gauge condition 
[47,36] in which we solve 
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(9) 



where At denotes a time step in numerical computation. 

During numerical simulations, violations of the Hamiltonian constraint and conservation of mass and angular 
momentum are monitored as code checks. Numerical results for several test calculations, including stability and 
collapse of nonrotating and rotating neutron stars, have been described in [35]. The axisymmetric code has been used 
for simulations of stellar core collapse to neutron stars and black holes, producing numerically convergent results [12]. 
The three-dimensional code has been used particularly for simulations of merger of binary neutron stars [41,36]. In 
[36], the details of the latest implementation are described, and we illustrate that accurate and convergent numerical 
results on the outcomes after the merger as well as on gravitational waveforms can be obtained with the present code. 



B. Equations of state 

A parametric equation of state is adopted following Miiller and his collaborators [6,10]. In this equation of state, 
one assumes that the pressure consists of the sum of polytropic and thermal parts as 

P = Pp + Pth. (10) 

The polytropic part is given by Pp = Kp(p)p^'^P^ where Kp and T are not constants but functions of p. This part 
corresponds to the cold (zero-temperature) part of the equation of state. In this paper, we follow [10] for the choice 
of Kp{p) and T{p): For the density smaller than the nuclear density which is defined as Pnuc = 2 x 10^* g/cm^, 
r = ri(=const) is set to be < 4/3, and for p > pnuc, T = r2(= const) > 2. Thus, 



KlP \ P< Pnuc, l^^Y) 



K2p'\ P>P, 



nuc? 



where Ki and K2 are constants. Since Pp should be continuous, the relation, K2 = Kip^-^^^^ , is required. Following 
[6,10], the value of Ki is fixed to be 5 x lO^'' in the cgs unit. With this choice, a realistic equation of state for p < Pnuc, 
in which the degenerate pressure of electrons is dominant, is approximated. Since the specific internal energy should 
be continuous at p = Pnuc, the polytropic specific internal energy ep is defined as 



£p = < 



Ki r 1 

P<Pn 



L r2-i^ + (ri-i)(r2-i) ' 



With this setting, a realistic equation of state for cold nuclear matter is mimicked for an appropriate choice of Fi and 
Fa. 

An advantage of the parametric equations of state is that we can investigate the dependence of the dynamics 
of stellar collapse on the equations of state systematically and very easily by changing the values of Fi and F2 
appropriately. A more realistic simulation with a realistic equation of state should be performed at the goal in this 
research field. However, the equations of state for p > pnuc are not still well-known. Also, because of complexity 
of the microphysical processes, in simiilations with such realistic equations of state, it is often not easy to extract 
essential physical properties of stellar core collapse such as key quantities that determine the maximum density in 
the collapse, the collapse time scale, the maximum value of T/W, the profiles of the density and angular velocity of 
formed protoneutron stars, T/W of formed protoneutron stars, nonaxisymmetric dynamical stabilities, and amplitude 
of gravitational waves. Simulations with the parametric equations of state are helpful to systematically answer these 
questions. 

In this paper, we choose (Fi,F2) = (1.3,2.5), (1.32, 2.25), and (1.28, 2.75). In Table I, we list the maximum mass 
and the corresponding density at the center for the three sets of Fi and F2 with pnuc = 2 x 10^^ g/cm^. In all 
three cases, the maximum ADM mass becomes about 1.6A'/q which is a reasonable value for neutron stars [48]. As a 
default, we set Fi = 1.3 and F2 = 2.5 in the following. In a previous Newtonian three-dimensional simulation [7], a 
different set as Fi = 1.28 and F2 = 2.5 is chosen, and the authors have found that only for such small value of Fi, 
nonaxisymmetric dynamical instabilities are induced. The choice of this set is acceptable in the Newtonian framework, 
but in general relativity it should not be adopted because with this choice, the maximum mass of a cold spherical 
neutron star becomes about I.3M0, which is too small for the maximum mass. Such choice should be excluded in 
general relativistic simulations. 
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M,{Mq) 


M{Mq) 


Pc(g/cm^) 


1.3 


2.5 


1.810 


1.600 


2.87el5 


1.32 


2.25 


1.754 


1.623 


2.39ol5 


1.28 


2.75 


1.869 


1.597 


3.18el5 


1.28 


2.5 


1.486 


1.298 


4.42el5 



TABLE I. Maximum baryon rest-mass, ADM mass, and the corresponding central density for spherical neutron stars of cold, 

parametric equations of state (11) for several choices of Fi and The first three sets of Fi and r2 arc adopted in the paper. 
The fourth set of the equation of state which is used in Refs. [6,7] is too soft and the maximum ADM mass for cold, spherical 
neutron stars is too small to be adopted as a realistic parameter set. 



The thermal part of the pressure Pth plays an important role in the case that shocks are generated. Pth is related 
to the thermal energy density eth = e — ep as 

Pth = (Fth - l)/0£th. (13) 

For simplicity, the value of Fth, which determines the strength of shocks, is chosen to be equal to Ti{^ 1.3). Our 
previous numerical work [12] showed that the results depend very weakly on the value of Fih as far as it is in the 
range between -^1.3 and 5/3. 

For the simulation, first, equilibrium rotating stars with F = 4/3 polytrope are given. Then, the simulations are 
started with equations of state (10). Since the value of the adiabatic index is slightly decreased from F = 4/3 to 
Fi(< 4/3), the collapse is triggered. The equilibrium states are computed adopting the polytropic equation of state 

P = Kop^^^ (14) 

where Kq is the adiabatic constant. In this paper, Kq is set to be 5 x 10^^, 7x 10^^, and 8 x 10^^ cm^/s^/g^/^. The latter 

two are adopted to increase the mass of the progenitor of stellar collapse: For the F = 4/3 polytrope, the mass (both 
the baryon rest-mass and the ADM mass) of the stars is approximately written as 4.555(iiro/G)'^/^ g, which depends 
very weakly on the rotational velocity profile [48] (cf. Table II). This implies that for if 14 = Kq/IO^^ cm'^/s^/g^/'^ = 5, 
7, and 8, the mass is about 1.5, 2.5, and 3Mq, respectively. Thus, for if 14 > 7, the total mass of the system is much 
larger than the maximum allowed mass of the cold spherical neutron stars chosen in this paper « 1.6Mq. 

For Kq ~ Kdeg « 5 X 10"'^'' cm^/s^/g^/"^ which is chosen in previous papers [6,10,12], a soft equation of state governed 
only by the electron degenerate pressure is approximated well [48,49]. On the other hand, the radiation pressure is 
also approximated by the F = 4/3 polytropic equation of state. Thus, by choosing Kq > ifdeg, we may consider that 
the pressure is composed of the sum of the electron degenerate pressure and the radiation pressure with the ratio 

Kdeg to ifrad = Kq - K^eg aS 

P = ifdegP^/'+^radp'/'. (15) 

In the simulation, Ki is related to ifdeg by Ki = ifdegPo^^ '"^ where we set po = 1 g/cm^. The specific internal 
energy is given by 

e = 3ifop'/^ (16) 

and the pressure at the initial stage is written as 

P = 3(Fi-l)ifop"/3, (17) 

implying that for Fi < 4/3, the pressure is depleted by (4 — 3Fi) = 4-16% for Fi = 1.32 1.28 at the initial stage. 
Namely, in this setting, with the smaller value of Fi, the pressure for a given value of p < pnuc becomes smaller, and 
also, the deletion factors of the pressure and the internal energy at the initial condition are larger. As shown in Sees. 
IV and V, the effect associated with the small change in Fi significantly modifies the dynamics of the collapse and 
the stability against nonaxisymmetric dynamical deformation. 

C. Wave extraction methods 

We extract gravitational waves using two methods. One is a gauge-invariant wave extraction method in which we 
perturbatively compute the Moncrief variables in a flat spacetime background [46] as we have used in our series of 
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papers (e.g., [45]). To compute them, first, we split 7jj into rjij + J^imdj^ spherical polar coordinates, where 

rjij is the flat metric and is given by 



/-Im. 

Si? 



/ H2lmyi'm hiimdgYi 



hllmdu>Yl 



r^KimYim + GimWim) T^GimXxr. 



* sin^ e{KimYim - GlmWlm) 

-CijndpYijn/ siri0 CijndeYijnSmO 

+ 1 * r'^DimXim/ sine -r'^DimWim sin 9 

* * -r'^DiniXim sin 9 



(18) 



Here, * denotes the symmetric components. The quantities i?2im,) hum, Kim, Gim, Cim, and Dim are functions of r 
and t, and are calculated by performing integrals over a two-sphere of a given coordinate radius [see [40] for details] . 
Yim is the spherical harmonic function, and Wim and Xim are 



Wlr. 



{def - cot Ode 



1 



sin^9 



Y, 



8$ — cot 9 



Yi, 



The gauge-invariant variables of even and odd parities are defined by 

R?mit, r)^^l ^ ( — + rdrD, 



i2{i-2y. 



2(1 + 2)1 fCi. 



{1-2)1 



where 



knm= Kim + l{l + i)Gim + 2rdrGim - 2 



hllm, 



k2lm = 



H2lm _ 

2 2dr 



r{Kim + l{l + l)Gim} 



(19) 

(20) 
(21) 

(22) 

(23) 



Using the gauge-invariant variables, the energy luminosity and the angular momentum flux of gravitational waves 
can be defined by 



dE 
Itt 



327r 



dJ r 
H ~ 32Tr 



Y,[\midtRfm)Rfm\ + \m{dtR?m)R?n 



l.m 



The total radiated energy and angular momentum are calculated by 



' , dJ 
dt — . 

dt 



(24) 
(25) 

(26) 



In this paper, we pay attention only to even-parity modes with I = 2 which are the dominant modes. 
To search for the characteristic frequencies of gravitational waves, the Fourier spectra are computed by 



RLif) = £ 



(27) 



where I denotes E and O. In the analysis, tf is chosen as the time at which the simulation is stopped. Before t < Vohs 
where robs denotes a radius at which gravitational waves are extracted, no waves propagate to robsj so that we choose 

Using the Fourier spectrum, the energy power spectrum is written as 



dE _ TT 2 



E fi\Rfm{fr + \R?mifr) (/>0), 



(28) 



;>2,m>0 
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where for m ^ 0, we define 



RL ^ yl\RLU)? + \Ri-mU)? (m > 0). 



(29) 



Note that in deriving Eq. (28), we use the relation |^fm(— /)| = l-Rfm(/)|- 

Computation of gravitational waves is also carried out in terms of a quadrupole formula which is described in 
[44,12]. As shown in [44], a kind of quadrupole formula can provide an approximate gravitational waveforms from 
oscillating compact stars. The quadrupole formula is in particular useful when the amplitude of gravitational waves 
is smaller than the numerical noise because in such case, it is difficult to extract gravitational waves from the metric 
in the wave zone. 

In quadrupole formulas, we compute gravitational waves from 



hij 



k r> I 



-PiiP''^ 

2 ^ 



r dt^ 



(30) 



where fij and Pi ' = Sij — niUj [ui = jr) denote a tracefree quadrupole moment and a projection tensor. 

In fully general relativistic and dynamical spacetimes, there is no unique definition for the quadrupole moment 
Following [44,12], we choose the formula as 



Then, using the continuity equation, we can compute the first time derivative as 



(31) 



(32) 



To compute , we carried out the finite differencing of the numerical result for 7^ . 

In this paper, we focus only on Z = 2 mass quadrupole modes. Then, the gravitational waveforms are described by 



■(1 + cos^ 6*) cos(2v?) + 4y(l + cos^ 6) sin(2(p) + /^^ - 



Ixx ~l~ ly 



'■yy 



cos 6 sin(2(^) + I^y cos 6 cos(2<^) 



in the quadrupole formula, and 



h+ = - 
r 



hx = 



^{R22+{1 + COS"" e) cos(2(^) + R22-{1 + cos^ e) sin(2(^)} + ^/ £^i?20 sin^ 9 
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— i?22+ COS 9 sin(2iy9) + i?22- COS 6 cos(2(p) 



in the gauge-invariant wave extraction where 

-R22± 



V2 



R20 = -Rfo''- 



(33) 
(34) 

(35) 
(36) 

(37) 
(38) 



For derivation of h+ and /ix, we assume that the wave part of the spatial metric in the wave zone is written as 

dl'^ = dr'^ + r^[(l + h+)d9'^ + sin^ 9(1 - h+)dip^ + 2si7i9hxd9dip], (39) 
and set I^z ~ lyz = and i?f ±1 = since we assume the reflection symmetry with respect to the equatorial plane. 



In the following, we present 



^+ ^xx ^yy ) 



2h 



"■yy 



(40) 
(41) 

(42) 
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in the quadrupole formula, and as the corresponding variables, 



R+ = ^ 


/ 5 „ 


(43) 


Rx = ^ 




(44) 


Ro = 


1 it^'°' 


(45) 



in the gauge- invariant wave extraction method. These provide the amplitude of a given mode measured by an observer 
located in the most optimistic direction. 



III. SETTING 



A. Initial conditions for axisymmetric simulation 

Rotating stellar cores in equilibrium with the F = 4/3 polytropic equation of state [Eq. (14)] are prepared for the 
initial conditions. Following [10,12], the maximum density is chosen as Pmax = 10^*^ g/cm^ irrespective of the velocity 
profile and the value of Kq. 

The velocity profiles of equilibrium rotating stellar cores are given according to a popular relation [50,51,10] 

= wjina - (46) 

where Q = v'^ denotes the angular velocity, fio is that along the rotational axis, and vud is a constant. In the 
Newtonian limit, the rotational profile is written as 

2 

^ = ^a-^f-^- (47) 

Thus, Wd indicates the steepness of differential rotation. Since the compactness of the initial data adopted in this 
paper is not so large with M/R ~ 10~^, where R denotes a stellar radius, that general relativistic effects are weak. 
As a result, the profile of the rotational angular velocity is approximately written by Eq. (47). In the following, we 
adopt rigidly rotating models in which oo, and differentially rotating models with A = vj^/Re = 0.25 and 

0.1, where i?e is the coordinate radius at the equatorial surface. The ratio of the angular velocity at the equatorial 
surface to Q.a is « 1/17 and 1/101 for A = 0.25 and 0.1, indicating that Eq. (47) is approximately satisfied. We pay 
particular attention to the case with high degrees of differential rotation in this paper, since in the collapse, a large 
value of [3 can be achieved only for such cases. Indeed, a study for presupernova evolution of rotating massive stars 
[38] indicates that the velocity profile of the iron core just before the onset of collapse may be differentially rotating. 

As introduced in Sec. I, the ratio (/? = T/W) of the rotational kinetic energy T to the gravitational potential 
energy W is often referred in the following. In general relativity, W is defined by [39] 

W ^ M^ + U + T ~ M. (48) 

Here, W is defined to be positive. For stable rotating stars in equilibrium with F = 4/3, is nearly equal to M. 
Thus, 

W^U + T, and P ^ (49) 

Even in the dynamical evolution, is the conserved quantity and M is approximately conserved in the case that 
luminosity of gravitational waves is small. Thus, if other components of the energy such as the kinetic energy associated 
with radial velocity are small, the above approximate relation for (3 in terms of U and T may be used even in the 
dynamical spacetime [see also discussion around Eq. (55)]. 

In Table II, several quantities for the models adopted in the present numerical computation are summarized. In the 
first column, we describe the name of each model. We refer to the models with {Kii,A) = (5,oo), (5, 0.1), (7, oo), 
(7, 0.25), (7, 0.1), (8,oo), and (8, 0.1) as M5a, M5c, M7a, M7b, M7c, M8a, and M8c, respectively. 

For the rigidly rotating case with A oo, the magnitude of the angular velocity has to be smaller than the Kepler 
angular velocity at the equatorial surface. This restricts the maximum value of /? to be smaller than 0.0089 for 
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TABLE II. Quantities for selected set of the initial conditions and the results of collapse are listed. Ko, M«, Po(= 27r/f2a), 
i?e, and PMax arc shown in units of cm'^/s^/g^/'', Mq, sec, km, and g/cm'^, respectively. Here, pMax and OMin arc the maximum 
and minimum achieved during the whole evolution. /3init and /Jmax denote the initial value of T/W and the maximum value 
of r/(r + U) achieved during the collapse. The baryon rest-mass M« is nearly equal to the ADM mass M for all the models. 
In the last column, the outcomes for Fi = 1.3, r2 = 2.5, and pnuc = 2 x 10^* g/cm^ are shown. Here, NS, O-A, and O-B 
denote that the outcomes are neutron star, oscillating star with the maximum density larger than pnuc, and oscillating star of 
subnuclear density, respectively. Note that for Ko = 8 x lO^** (cgs) and A — »■ oo, any star collapses to a black hole and that for 
Ko < 6 X 10^** (cgs), any star does not collapse to a black hole. 



all the values of Kq. This implies that the angular velocity for models M5al, M7al, and M8al is approximately 
maximum among the rigicily rotating cases for a given value of Kq. The final outcome of M8al is a black hole. This 
implies that any star with {Ki4,A) = (8, cx)) collapses to a black hole because the mass is too large and the angular 
momentum is too small to halt the collapse. On the other hand, for if 14 < 6, any star does not collapse to a black 
hole since the mass is not large enough. The detail on the criterion of the formation of black hole is also described in 
the companion paper [52]. 

B. Method of axisymmetric simulation 

During the collapse, the maximum density increases from 10^*^ g/cm^ to ~ 10^^ g/cm^ in the neutron star formation 

and to ^ 10^^ g/cm^ in the black hole formation. This implies that the characteristic length scale of the system varies 
by a factor of ^ 100. In the early phase of the collapse which proceeds in a nearly homologous manner, we may follow 
the collapse with a relatively small number of grid points by moving the outer boundary inward while decreasing 
th(^ grid spacing, without increasing the grid number by a large factor. As the collapse proceeds, the central region 
shrinks more rapidly than the outer region does and, hence, a better grid resolution is necessary to accurately follow 
such a rapid collapse in the central region. On the other hand, the location of the outer boundaries should not be 
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changed by a large factor to avoid discarding the matter in the outer envelopes. 

To compute such a collapse acciiratcly while saving the CPU time efficiently, a regridding technique as described 
and used in [53,12] is adopted. The regridding is carried out whenever the characteristic radius of the collapsing star 
decreases by a factor of a few. At each regridding, the grid spacing is decreased by a factor of 2. All the quantities in 
the new grid are calculated using the cubic interpolation. To avoid discarding the matter in the outer region, we also 
increase the grid number at the regridding, keeping a rule that the discarded baryon rest-mass has to be less than 1% 
of the total. 

Specifically, N and L in the present work are chosen using a relativistic gravitational potential defined by $c = 
1 — ttc ($c > 0), which is ^ 0.01 at t = 0. Here, ac denotes the central value of the lapse function. Since <I>c is 
approximately proportional to M/R, can be used as a measure of the characteristic length scale for the regridding. 
Typically, the value of N is chosen monitoring the magnitude of $c in the following manner; for $c < 0.04, we set 
N = 420; for 0.04 < $c < 0.1, we set N = 700; for 0.1 < *c < 0.2, we set N = 1200; and for 0.2 < $c, we set 
N — 1800, and keep this number until the termination of the simulations. Note that at t = 0, the equatorial radius 
is covered by 400 grid points in this case. With this setting, the total discarded fraction of the baryon rest-mass 
which is located outside new regridded domains is ^ 1%. The grid spacing in iV = 1800 is ~ 0.6 km for differentially 
rotating initial condition and ~ 0.6-0.8 km for rigidly rotating cases. A previous work [12] illustrates that with these 
grid resolutions a convergent result is obtained for most cases. 

Nevertheless, we should be very careful in judging black hole formation since the criterion for the black hole 
formation near a threshold depends sensitively on the strength of shocks that are formed when the density around the 
central region exceeds pnuc- The shocks in numerical simulations in general become stronger with improving the grid 
resolutions. This implies that a black hole may be spuriously formed in a coarse grid resolution in which the strength 
of the shocks is underestimated. To avoid such misjudging, in the case that a black hole is likely to be formed, we 
perform simulations with a finer grid resolution as follows; for $c < 0.04, we set N = 620; for 0.04 < <I>c < 0.1, we 
set N = 1020; for 0.1 < <I>c < 0.2, we set N = 1700; and for 0.2 < we set N = 2500. Note that at t = 0, the 
equatorial radius is covered by 600 grid points in this case. If we find a convergent result on the black hole formation 
in both resolutions, we judge that the black hole is formed. 

Simulations for each model with the typical grid resolution are performed for 40,000-50,000 time steps. The required 
CPU time for one model is about 20 hours using 4 processors of FACOM VPP5000 at the data processing center 
of National Astronomical Observatory of Japan, and about 10 hours using 8 processors of NEC SX6 at the data 
processing center of ISAS in JAXA. 

C. Method of three-dimensional simulation 

Since computational resources are restricted and we cannot take the grid number per one direction as large as that 
in the axisymmetric case, it is not a good idea to perform three-dimensional simulations from the initial conditions 
with pmax = 10^^° g/cm^. To save computational time, we always follow the early stage of the collapse using the 
axisymmetric code. After the collapse proceeds sufficiently, we change to the three-dimensional code. In preparing 
the initial conditions of three-dimensional computations, numerical results of the axisymmetric simulations are used. In 
the early stage of collapse at which the value of /? of the collapsing star is not still very large 0.2), nonaxisymmetric 
dynamical instabilities will not be induced. For highly differentially rotating cases, nonaxisymmetric instabilities could 
be induced even with a low value of /3 [28]. However in such cases, the growth time of the nonaxisymmetric instabilities 
would be much longer than the collapse time scale [28]. Therefore, the method that we adopt is appropriate. 

Specifically, the initial condition for the three-dimensional simulations is prepared when the central value of the 
lapse function becomes ac — 0.8 in the axisymmetric simulations. (For some case in which the minimum value of ac 
is slightly larger than 0.8, we choose the value as 0.85.) Since our major purpose in the three-dimensional simulations 
is to investigate the nonaxisymmetric dynamical stability of the collapsing star, we add a nonaxisymmetric density 
perturbation to the axisymmetric state. Associated with this change, the metric should be also perturbed, but we do 
not know how to do. For this reason, we adopt a very simple method for setting the initial conditions as follows. 

First, we note that for ac > 0.8, the magnitude of jij — Sij is very small (<C 0.01) for all the components, and hence, 
the spatial hypersurfacc is approximately conformally flat. Also, the trace of the extrinsic ciirvature is nearly equal 
to zero because of our choice of the slicing condition. Thus, in setting the initial conditions of the three-dimensional 
simulations, we assume that the three-hypersurface is conformally flat and K/J' = for simplicity. Then, we determine 
the conformal factor and the trace- free extrinsic curvature using the constraint equations. In this case, for a solution 
of the constraint equations, we only need to extract p*, e, and Ui from the numerical results of the axisymmetric 
simulations in the following method. 

In the first step, we solve the momentum constraint equation using the York's procedure [54]. Setting 
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Aij = Aiji,^ = diWj + djWi - ^6ijdkWk, (50) 

where Wk is a three-vector, the momentum constraint is written as 

Aflat Wi + ^didjWj = Snp^Ui, (51) 

where Anat is the flat Laplacian. Since and m, are given, this equation is solved by a standard procedure (e.g., 
[13]) to give the tracefree part of the extrinsic curvature. 

In the next step, the Hamiltonian constraint equation is solved. In the conformally flat spatial hypersurface, the 
equation for the conformal factor ip is written as 

Aflat^ = -2TTp,etp-^ - -—A,jA'^ . (52) 

o 

Since and e were given, and Aij was already computed in the first step, ip is also computed by a standard procedure 
often used in the initial value problem. 

The simulations were performed using a fixed uniform grid and assuming reflection symmetry with respect to the 
equatorial plane. The typical grid size is (27V + 1,2A^ + 1, A'^ + 1) for {x,y,z), and we adopt N = 156, 188, and 
220. The grid covers the region —L < x < L, —L < y < L, and < z < L where L is the location of the outer 
boundaries along each axis. For a given model, we take the identical value of L irrespective of the value of N. The grid 
spacing Ax = L/N is chosen to be larger than that adopted in the corresponding axisymmetric simulation because of 
restricted computational resources for the three-dimensional case. In the case of N = 156, we choose the grid spacing 
twice as large as that of the corresponding axisymmetric simulations. The typical computation is performed with 
A'^ = 188, and to check the convergence, the value of N is varied. For models in which a bar-mode instability sets in, 
simulations are performed with N = 220. 

The value of L is much smaller than that of the axisymmetric simulation. This implies that we discard the matter 
located in the outer region of the collapsing core. Specifically, we discarded the matter outside a sphere of radius « tq 
by the rule 

p*(new) = /3* (axisymmetric) x _^ ^ , (53) 

where ro « 0.95L and 5r = 2Ax. In this method, the fraction of the discarded baryon mass located for r > ro is 
about 10 -20% (compare Tables II and IV). 

In this paper, we focus primarily on the dynamical stability against m — 2 bar-mode deformation, since it is 
expected to be the fastest growing mode. Specifically, we superimpose a density perturbation in the form 

= p*(new) (^1 + 0.4 ^ "/ ^ . (54) 

To check that the bar-mode perturbation grows for dynamically unstable models even when the initial configuration 
is nearly axisymmetric, we also performed simulations without adding nonaxisymmetric perturbation besides random 
numerical noises for selected unstable models. We found that in such cases, the bar-mode perturbation indeed grows 
although it takes more computational time to follow the growth. 

In the case that the equations of state are very soft, the degree of the differential rotation is very high, and the 
value of f3 is large enough (> 0.14) for a collapsed star, m = 1 modes may grow faster than the m = 2 mode [55,56]. 
In the formation of neutron stars in which pmax > Pnuc, the equation of state is stiff, and hence, the m = 1 mode 
may not be very important. On the other hand, in the formation of oscillating stars, equations of state can be soft 
for pmax < Pnuc- Howcver, the value of (3 in such phase of subnuclear density are not very large. Thus, it is expected 
that even if the m = 1 mode becomes unstable, the perturbation may not grow as significantly as found in [55,56]. 
Hence, we do not pay particular attention to this mode in this paper. Since nonaxisymmetric numerical noises are 
randomly included at t = 0, in some models, the m = 1 mode grows as found in Sec. V. However, the amplitude of 
the perturbation is indeed not as large as that for m = 2. 

Since we assume the conformal flatness in spite of the fact that the conformal three- metric is slightly different from 
zero in reality, a small systematic error is introduced in setting the initial data. Moreover, we discard the matter 
located in the outer region of the collapsing core according to Eq. (53). This could also introduce a systematic error. 
To confirm that the magnitude of such error induced by these approximate treatments is small, we compare the results 
in the three-dimensional simulations with those in the axisymmetric ones. We have found that the results agree well 
each other and the systematic error is not very large. This will be illustrated in Sec. V (cf. Fig. 13). 
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Simulations for each model with the grid size (441, 441, 221) (TV = 220) were performed for about 15,000 time 
steps. The required CPU time for computing one model is about 30 hours using 32 processors of FACOM VPP 5000 
at the data processing center of National Astronomical Observatory of Japan. 



IV. NUMERICAL RESULTS OF AXISYMMETRIC SIMULATIONS 




A. Outcomes 



In the last column of Table II, we summarize the outcomes of stellar core collapse in the axisymmetric simulations 

for Fi = 1.3 and F2 = 2.5. They are divided into three types: Black hole, neutron star, and oscillating star for which 
the maximum density inside the star is not always larger than Pnuc- For given values of Ko{> 7 x 10^'* cgs) and A, a 
black hole is formed when the initial value of P (hereafter Anit) is smaller than critical values that depend on A. As 
described in Sec. Ill A, f3 in the collapse is defined by 

In the dynamical spacetime with M* RiMforF = 4/3,W would be approximately written as 

W^U + T + Tother, (56) 

where Tother denotes a partial kinetic energy obtained by subtracting the rotational kinetic energy from the total. 
Thus, T/W should be approximated hy T/{U + T + Tocher), but we do not know how to appropriately define Tother- 
Fortunately, it would be much smaller than T at the initial state, at the maximum compression at which the spin 
of the collapsing star becomes maximum, and in a late phase when the outcome relaxes to a quasistationary state. 
This implies that using the definition of (55), the maximum value and a final relaxed value of (3 will be computed 
approximately. In other phases, (3 computed by Eq. (55) gives an overestimated value. 

In Fig. 1, we show the evolution of the central value of the lapse function (etc), the maximum value of the density 
(Pmax), and (3 for models M7cl, M7c2, M7c3, M7c5, and M7c6. In the following, we denote the maximum density 
and minimum value of the lapse achieved in the whole evolution as pMax and ajviin, respectively. On the other hand, 
the maximum density at a given time is denoted by Pmax- 
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FIG. 2. (a) QfMin and pMax and (b) /3max for various values of /3init. In both panels, the solid triangle, the solid circles, open 
squares, open triangles, and crosses denote models M5c, M7a, M7b, M7c, and M8c, respectively. The dotted line in panel (b) 
denotes /3max ~ 0.27. 



Figure 1 shows that for most cases, the value of pMax becomes larger than pnuc- However, with the increase of /?init, 
it decreases significantly. Also, for several cases, Pmax drops below p^uc soon after it reaches the maximum. Such 
oscillating stars for which the values of Pmax oscillate between pi(> Pnuc) and P2(< Pnuc) are referred to type 0-A 
in Tabic IL On the other hand, if /3init is not very large and neither is the maximum value of f3 (hereafter /3max), a 
neutron star or a black hole is formed. Here, formation of a neutron star implies that Pmax achieved after the stellar 
collapse is always larger than pnuc- Formation of a black hole implies that we confirm the formation of apparent 
horizon. 

In Figs. 2(a) and (b), we show OfMin, PMax, and /3max for various values of /Sinit- For /Jinit ^ 0.02 with Kq = 7 x lO^** 
cgs, PMax is smaller than pnuc, and the resulting star is quasiradially oscillating with the subnuclear density. Such 
stars are referred to as 0-B in Table II. 

Figure 1 and Table II show that initial high degrees of differential rotation with A ~ 0.25 and 0.1 have an effect for 
preventing black hole formation. (Compare the parameters among the models with Kq — 7 x 10^^ cgs.) For the rigidly 
rotating cases, the stars with /Jjnit ^ 0.005 collapse to a black hole. On the other hand, the stars with /3init ~ 0.003 
do not collapse to a black hole but form a neutron star for A = 0.1 and 0.25. This is simply because the stars with 
such high degrees of differential rotation have a large centrifugal force near the rotational axis, and hence, even in 
the case that the global value /3init is small, the effective local value of the centrifugal force would be large enough to 
prevent cores from collapsing to a black hole. 



B. Evolution of /3 for Ti = 1.3 and r2 = 2.5 

Figures 1 and 2 indicate that with the increase of /?init, ctMin (pMax) increases (decreases). The value of /?max is 
larger for the larger value of /3init as far as /3init ^ 0.02. However, the amplification factor /?max/Anit is smaller for the 
larger value of (Sjuit- This is because in the collapse with the large values of /3init, strong centrifugal force prevents the 
maximum value of compactness (or maximum density or maximum value of the gravitational potential) from being 
increased by a large factor. Spin can be increased by a larger factor for a star which gains a larger compactness. 
Therefore, for stars of approximately identical mass, the amplification factor /3max/Anit should be smaller for the 
larger value of /3init • 

The typical value of /3max/Anit is 10-20 for rigidly rotating case and for diflFerentially rotating cases with /3init ^ 0.015. 
Naive estimation predicts that T oc J^/(Mi?^) and W cx M'^/R, and hence, /? (x J'^/{M^R). Thus, /? seems to be 
proportional to the inverse of the stellar radius in the condition that the mass and the angular momentum of the 
system are conserved. Since the characteristic stellar radius changes by a factor of ^ 100 during the collapse, we 
may predict that [3 also increases by two orders of magnitude. However, this does not occur. The reason for the 
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(c) (d) 
FIG. 3. (a) and (b): The density contour curves for p for model M7cl at i = 116.0 and 118.3 msec. The solid contour curves 
are drawn for p/ p^nax^ = e~°'^^ for = 1, 2, 3, • • • , 20 where Pmax denotes the maximum of p at the given times, which are found 
from Fig. 1(a). Vectors indicate the local velocity field and the scale is shown in the upper right-hand corner, (c) 

Density profiles along the x and z axes aX t = 116.0 (solid curve) and 118.3 msec (dashed curve), (d) The same as (c) but for 
angular velocity along the radial coordinate in the equatorial plane. 



rigidly rotating case is that although the core radius decreases to ~ 10 km in the collapse, the outer region of the star 
which possesses a large fraction of the angular momentum does not collapse to such a small radius due to the strong 
centrifugal force. The reason for the highly differentially rotating cases with a high value of /3init ^ 0.015 is that the 
centrifugal force near the rotational axis is so strong that the collapse is halted before the stellar radius becomes ~ 10 
km. For low values of /3init, ^ 0.01, the rotational velocity in the outer region is small, and also, the centrifugal force 
in the central region is not as strong as that for /3i„it > 0.015. As a result, the stellar components that enclose a large 
fraction of the angular momentum can collapse to small radii, and hence, /3 can increase by a factor of > 30. 

For given values of Kq and /9init, the value of /3niax is larger for higher degrees of differential rotation. This suggests 
that stellar cores with a higher degree of differential rotation may be more subject to nonaxisymmetric dynamical 
instabilities. For the identical value of Anit, /3max is larger for higher-mass stellar cores with a larger value of Kq in 
the case of highly differentially rotating cores. (Compare, e.g., the solid triangle, the open triangle, and the cross in 
Fig. 2(b); more specifically, compare the results for models M5cl, M7c2, and M8c3, for which the values of /3init are 
approximately the same as 0.0177, but /3niax is larger for larger mass.) The reason for this behavior is that the stars 
of higher mass can reach a more compact state during the collapse, and as a result, their spins can be increased by a 
larger factor and so can be (3. On the other hand, for rigidly rotating cases, this feature is not very outstanding. 

An interesting point is that the value of /3max has a maximum around /3i„it ^ 0.018 for A = 0.1 and 0.25: For 
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FIG. 4. The same as Fig. 3(a)-(d) but for model M7c3 a.t t = 114.9 and 120.4 msec. The solid and dashed curves in panels 
(c) and (d) are drawn for the corresponding time slices, respectively. 



Anit ~ 0.018, /3max IS an increase function of Anit) reflecting the initial magnitude of the spin. However, for /3init ^ 
0.018, /3max is a decrease function. The reason is that the centrifugal force of the rotating stars with /3init ^ 0.018 is 
so strong that the collapse is halted before the stellar core becomes compact enough. This feature is also reflected in 
Fig. 2(a) from which we find that the value of aMin is an increase function of /3init- 

As reviewed in Sec. I, nonaxisymmetric dynamical instabilities of rotating stars in equilibrium set in when the value 
of /3 becomes larger than ~ 0.27. If we assume that the collapsing stars with /3max ^ 0.27 are dynamically unstable. 
Fig. 2(b) suggests that the conditions for the onset of the instabilities will be the following: (i) the progenitor of 
the collapse should be highly differentially rotating with A < 0.25; (ii) the progenitor has to be moderately rapidly 
rotating with 0.01 < /3init ^ 0.02; (iii) the progenitor star has to be massive enough. 

However, it should be kept in mind that the condition /3max > 0.27 is satisfied only for a few msec during the stellar 
collapse. This indicates that if the growth time scale of nonaxisymmetric instabilities is not as short as a few msec, 
the system may remain nearly axisymmetric. Thus, the condition /Smax ^ 0.27 does not have to be the criterion for 
the onset of nonaxisymmetric dynamical instabilities in dynamical systems. The examples are shown in Sec. V. 



C. Profiles of density and angular velocity for Fi = 1.3 and F2 = 2.5 

In [21], Tohline and Hachisu illustrated that the stars with toroidal density profiles are dynamically unstable against 
a bar-mode perturbation even if /? is much smaller than 0.27. Also in [28], we indicated that not only the value of /? 
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FIG. 5. The same as Fig. 3(a)-(d) but for model M7c6 a.t t = 113.0 and 117.0 msec. The solid and dashed curves in panels 
(c) and (d) are drawn for the corresponding time slices, respectively. For comparison, the results for model M7a2 at t = 121.3 
msec are shown (dotted curves). 



but also the degree of differential rotation is a key parameter for determining nonaxisymmetric stability of rotating 
stars. Thus, here, we focus on the profiles of the density and the rotational angular velocity of the outcomes in the 
stellar collapse. 

In Figs. 3-5, we display the snapshots of the density contour curves, velocity vectors in the x-z plane, density profiles 
along X and z axes, and rotational angular velocity profiles as a function of the radial coordinate in the equatorial plane 
at the time slices that the maximum completion is achieved and the system relaxes to an approximately quasistationary 
state for models MTcl, M7c3, and M7c6. These figures clarify how the outcomes are changed with the varying /?init 
for (approximately) identical values of A and M. Panels (a), (b), and (c) show that for the larger values of Anit) the 
shape of the outcome is more torus-like. Numerical studies for nonaxisymmetric dynamical instabilities in rapidly 
rotating stars in equilibrium have illustrated that torus-like stars arc often unstable [21,22,28]. This indicates that the 
models such as M7cl-M7c3 in which torus-like structures are formed are candidates for the onset of nonaxisymmetric 
dynamical instabilities even when the value of f3 is smaller than ~ 0.27. 

The panels (d) in Figs. 3-5 show that all the outcomes of the collapse arc differentially rotating. The degree of 
the differential rotation is very large for the cylindrical radius tu <^ 10 km as oc zu^^ with 5 ^ 1.9 2.0, reflecting 
the initial profile. In the inner region of tu ^ 10 km, the rotational angular velocity does not change as steeply as 
that for w ^ 10 km. This also seems to reflect the initial rotational velocity profile for which O is nearly constant 
for zu ^ zud- However, except for the very inner region, the star is totally diflferentially rotating, in particular in the 
outer region, for any models of A = 0(0. 1). The results indicate that the initial rotational velocity profile is reflected 
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0.28 


O-A 


M5clC 


1.28 


2.75 


0.80 


3.8cl4 


0.33 


0-A 


M5c2A 


1.3 


2.5 


0.74 


5.0el4 


0.28 


0-A-+NS 


M5c2C 


1.28 


2.75 


0.77 


5.2el4 


0.31 


0-A^NS 


M7b3A 


1.3 


2.5 


0.65 


4.8el4 


0.29 


NS 


M7b3B 


1.32 


2.25 


0.71 


2.8el4 


0.23 


O-A 


M7b3C 


1.28 


2.75 


0.70 


5.3el4 


0.29 


NS 


M7c2A 


1.3 


2.5 


0.61 


6.0el4 


0.33 


O-A 


M7c2B 


1.32 


2.25 


0.75 


5.6el4 


0.21 


O-A 


M7c2C 


1.28 


2.75 


0.67 


5.6el4 


0.36 


O-A^NS 


M7c3A 


1.3 


2.5 


0.51 


9.2el4 


0.30 


NS 


M7c3B 


1.32 


2.25 


0.38 


1.2el5 


0.26 


O-A 


M7c3C 


1.28 


2.75 


0.61 


7.6el4 


0.33 


NS 


M7c4A 


1.3 


2.5 


0.46 


1.2el5 


0.27 


NS 


M7c4B 


1.32 


2.25 








BH 


M7c4C 


1.28 


2.75 


0.59 


8.7el4 


0.30 


NS 


M8c2A 


1.3 


2.5 


0.59 


5.4el4 


0.34 


O-A 


M8c2B 


1.32 


2.25 


0.86 


2.2el4 


0.16 


O-A 


M8c2C 


1.28 


2.75 


0.64 


5.2el4 


0.37 


O-A 


M8c4A 


1.3 


2.5 


0.29 


1.5el5 


0.30 


NS 


M8c4C 


1.28 


2.75 


0.51 


9.4el4 


0.34 


NS 



TABLE III. Numerical results with different values of Fi and r2 for selected models. 



in the outcomt;. Thus, if the progenitor of the coUapse is highly differentially rotating, the outconies will be always 
so and, as a result, be candidates for the bar- mode dynamical instabilities [28]. 

In Figs. 5(c) and 5(d), we display together the profiles of the density and the rotational angular velocity for 
model M7a2 (dotted curves) at i = 121.3 msec at which the outcome has already relaxed to a quasistationary state. 
Figure 5(c) shows that the outcome is a spheroid not a torus-like object (i.e., the central density is highest). This 
is a characteristic property in the collapse with rigidly rotating initial conditions [52]. Figure 5(d) shows that the 
rotational angular velocity in the inner region is approximately flat and thus the high-density part of the protoneutron 
star is approximately rigidly rotating. The outer region of lu > 10 km, on the other hand, is differentially rotating, 
but the rotational angular velocity falls off in proportional to with 5 ^ 1.4 1.5: i.e., the profile is approximately 
that of Kepler's law, and hence, the degree of differential rotation in this case is smaller than that for differentially 
rotating initial conditions. More details about the outcomes in the rigidly rotating initial conditions are found in [52]. 

From the density contour curves, it is found that for all the differentially rotating models, the column density 
integrated along the rotational axis is much smaller than that along the equatorial plane after shocks sweep the 
matter. Namely, a funnel is formed around the rotational axis even in the absence of a black hole. This is due to the 
facts that the total mass around the rotational axis is initially small because of a high degree of differential rotation 
for the initial condition and that the formed shocks are strongest around the rotational axis. A current popular model 
for the central engine of gamma-ray bursts is the so-called collapsar model [57]. To escape the baryon- loading problem 
for the fire-ball model [58], it is often required to form a funnel in the collapsar models. In their scenario, a rapidly 
rotating black hole is formed, and subsequently, a jet emitted along the rotational axis of the black hole ejects the 
matter. The present results suggest that a high degree of differential rotation for the progenitor of the stellar collapse 
helps making a funnel without relying on the formation of a rapidly rotating black hole and subsequent jets. 

D. Dependence on equations of state 

To clarify the dependence of the outcomes on the equations of state, we performed simulations varying Fi and F2 
as hsted in Table III for models M5cl, M5c2, M7b3, M7c2, M7c3, M7c4, M8c2, and M8c4. Here, we focus only on 
highly differentially rotating cases with A = 0.1 and 0.25. The details for the cases of rigid rotation and moderate 
degrees of differential rotation with Ar^\ are presented in [52]. As listed in Table I, we choose three sets of (ri,r2) 
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as (i) (1.3, 2.5), (ii) (1.32, 2.25), and (iii) (1.28, 2.75). In the following, we will refer to the models with (i), (ii), and 
(iii) using the labels A, B, and C, e.g., as M5clA, M5clB, and M5clC. 
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FIG. 6. Evolution of (a) ac and (b) /3 for model M7c3 and of (c) ac and (d) /? for model M5c2. Panels (i), (ii), and (iii) are 
results for (ri,r2) = (1.3,2.5), (1.32, 2.25), and (1.28, 2.75), respectively. The dotted horizontal lines in panels (b) and (d) 
denote /3 = 0.27. 



In Fig. 6, we compare the evolutions of the central value of a and /3 for models M7c3A, M7c3B, and M7c3C 

(Figs. 6(a) and (b)) and for M5c2A and M5c2C (Figs. 6(c) and (d)) as representative illustrations. In the previous 
section, we found that models M7c3 and M5c2 are possible candidates for the onset of nonaxisymmetric dynamical 
instabilities of /3max > 0.27. Among these models of different equations of state, /3max for case (iii) is largest. On the 
other hand, /3max for case (ii) is much smaller than those in other two cases. This seems to be due to the fact that 
for smaller values of Fi, the depletion factor of the internal energy and the pressure in an early stage of collapse in 
which p <C Pnuc is larger. As a consequence, the internal energy U is decreased to contribute to the increase of /3, and 
furthermore, the matter around the rotational axis which possesses large values of the specific angular momentum 
collapses to a more compact state, for which a spin-up is enhanced effectively. On the other hand, the values of ac 
(compactness) for models M7c3C and M5c2C are larger (smaller) than that for M7c3A and M5c2A, respectively. This 
may be partly due to the fact that r2 for case (iii) is larger than that for (i), but mainly to the fact that the fraction 
of the matter which simultaneously collapses is smaller for case (iii) than for (i). This implies that to achieve a large 
value of (3, it is not necessary for the whole system to become compact. Rather, essentially needed is to accumulate 
the matter with large values of the specific angular momentum in the central region. This point is reconfirmed from 
the results for M7c3B. In this case, the value of /3max is much smaller than those for other two cases, although the 
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value of ttc is smallest among three eases. This is due to the fact that in this case, a large fraction of the matter 
collapses nearly simultaneously independent of the magnitude of the specific angular momentum. 

Table III also shows that the largest value of /3,nax is achieved for case (iii) for all the initial conditions. This 
indicates that to achieve a large value of (5, the depletion of the internal energy and the pressure in the early stage of 
the collapse, which in reality will be achieved by partial photo-dissociation of the iron to lighter elements and by the 
electron capture [48,49], should bo sufficiently large to accelerate the collapse of the central region. 

As indicated in Fig. 6, the outcomes for models M7c3 and M5c2 depend sensitively on the equations of state. For 
the small values of Fi (cases (i) and (iii)), an oscillating protoneutron star is formed eventually. The amplitude of the 
oscillation is smaller and the period is shorter for the smaller value of Fi (case (iii)). As a result, the protoneutron 
star relaxes to a quasistationary state more quickly. For a long period of the oscillation, the duration of the phase, 
in which /? and rotational angular velocity remain small, becomes long. This also suggests that for the onset of 
nonaxisymmetric dynamical instabilities, the smaller value of Fi may be preferable. 
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FIG. 7. The same as Fig. 3(a)-(d) but for model M7c3C at i = 85.1 and 89.0 msec. The solid and dashed curves in panels 
(c) and (d) are drawn for the corresponding time slices, respectively. In panels (c) and (d), we plot the results for model M7c3A 
at t = 120.4 msec (dotted curves) for comparison. 



In Fig. 7, we display the snapshots of the density contour curves and the velocity vectors at the time slices that 
the maximum compression is achieved and the system relaxes to an approximately quasistationary state for model 
M7c3C. The density profiles and the angular velocity in the equatorial plane at the corresponding time steps are 
shown in panels (c) and (d) together with results for model M7e3A (displayed in Fig. 4). It is foimd that the shape 
of the outcome for M7c3C is more torus-like than that for M7c3A. In addition, the degree of differential rotation for 
M7c3C is slightly higher than that for M7c3A. These facts indicate that the outcome for M7c3C is likely to be more 



19 



subject to nonaxisymmetric dynamical instabilities. These properties depend weakly on the value of Kq (i.e., mass 
of the progenitor). Indeed, the outcome of M5c2C is more torus-like and the degree of differential rotation in the 
central region for M5c2C is higher than those for M5c2A. As illustrated in Sec. V, thus, the value of Fi is one of key 
parameters for determining the onset of nonaxisymmetric dynamical instabilities. 




X (km) t (msec) 

(c) (d) 
FIG. 8. (a) The density contour curves for p for model M7b3A at t = 123.1 msec. The contour curves and vectors arc drawn 
in the same method as in Fig. 3. (b) The same as panel (a) but for model M7b3C at t = 91.3 msec, (c) Angular velocity along 
the radial coordinate in the equatorial plane for models M7b3A at t = 123.1 msec (solid curve) and M7b3C at t = 91.3 msec 
(dashed curve), (d) pmax as a function of time for models M7b3A, M7b3B, and M7b3C. 

We note that the properties pointed above is found only for A = 0.1. For A = 0.25, the value of /3max for two 
equations of state (i) and (iii) are not very different. Also, the density profile and the angular velocity profile of the 
formed neutron stars are similar (see Fig. 8). This indicates that for large values of A, the larger depletion factor 
of the internal energy (smaller value of Fi) in the early stage of the collapse does not play an important role for 
accumulating the matter of large specific angular momentum in the central region. This result suggests that the 
stability property against nonaxisymmetric deformation will not depend on the choice of Fi and F2 for A = 0.25 as 
strongly as for A = 0.1 as long as Fi < 1.3. On the other hand, for Fi = 1.32, the outcomes are completely different 
from those for other two cases as in the case of ^ = 0.1. 
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TABLE IV. Parameters and numerical results for three-dimensional simulations. Ko, M^(M), pmax, and L are listed in 
units of cgs, Mq , g/cm^ , and km, respectively. The ADM mass M is still nearly equal to the baryon rest-mass M* . In the last 
column, the stability against the bar-mode is shown. 



E. Candidates of noncLxisymmetric dynamical instabilities 

As reviewed in Sec. I, nonaxisymmetric dynamical instabilities of rotating stars in isolated equilibrium may set 
in when the value of /3 becomes larger than ~ 0.27 or when the degree of differential rotation is sufficiently high. 
It is found that to achieve /3max ^ 0.27, the following conditions are necessary; (i) the progenitor of the collapse 
should be highly differentially rotating with A < 0.25; (ii) the progenitor has to be moderately rapidly rotating with 
0.01 < /3init ^ 0.02; (iii) the progenitor should be massive enough to make a compact core for which an efficient 
spin-up is possible. 

As indicated in [21,28], even in the case of /3 < 0.27, nonaxisymmetric dynamical instabilities may set in if the 
degree of differential rotation is sufficiently large. To achieve such a state, the conditions (i) and (ii) are also necessary. 

In addition, the following condition is required: (iv) the depletion factor of the internal energy and the pressure in 
an early stage of collapse during which p <^ p^uc is large enough to induce a significant collapse in the central region 
for making a torus-like structure and a steep profile of rotational angular velocity. In the next section, we present 
numerical results of the three-dimensional simulations and illustrate that the condition (iv) plays an important role 
for the onset of nonaxisymmetric dynamical instabilities. 

V. RESULTS IN THREE-DIMENSIONAL SIMULATIONS 
A. Features of nonaxisymmetric dynamical instabilities 

As described in Sec. IV, there are several candidate models for which nonaxisymmetric dynamical instabilities may 
set in. We performed the three-dimensional simulations focusing on the candidates in which /3max ^ 0.27 and (3 ~ 0.2 
after bounce. In this paper, however, we do not pay attention to the models in which oscillating stars with the period 
^ 10 msec arc formed since the simulations for such models take too much computational time. 

In the simulations, we initially superimposed a nonaxisymmetric bar- mode density perturbation as defined in Eq. 
(54). Specifically, we picked up models M5cl, M5c2, M7b3, M7c2, M7c3, and M7c4 with (ri,r2) = (1.3,2.5) and 
(1.28, 2.75) (referred to them, e.g., as M5clA and M5clC) as listed in Table IV. Since the matter in the outer region 
is discarded in preparing the initial conditions for the three-dimensional simulations according to Eq. (53), the mass 
and the angular momentum are smaller than those in the corresponding axisymmetric simulations by 10-20%. As a 
consequence, the numerical results deviate slightly from those obtained by the axisymmetric simulations even in the 
case that nonaxisymmetric deformation is small. However, qualitative differences between two results are not found 
and also the quantitative disagreement is small (see below). 

In this paper, the dynamical stability against bar-mode deformation is analyzed using a distortion parameter defined 
by 

rj^irjl + rjir/^ (57) 

where 
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FIG. 9. Evolution of pmax, OLc, and r] in the three-dimensional simulations (a) for models M5clC (solid curves), M5c2A (dashed 
curves), and M5c2C (dotted curves), (b) for models M7c2A (dotted curves), M7c2C (solid curves), M7c3A (dotted-dashed 
curves), and M7c3C (dashed curves), and (c) for models M7b3A (dotted curves) and M7b3C (solid curves). Here, t2d denotes 
the time at which we change to the three-dimensional code; tiA =85.4, 115.1, 85.1, 114.9, 85.1, 121.2, 84.7, 118.6, and 86.7 
msec for M5clC, M5c2A, M5c2C, M7c2A, M7c2C, M7c3A, M7c3C, M7b3A, and M7b3C, respectively. 



Qxx Qyy '^Qxu 

^ ^yy ^xx ^ ^yy 

and 

Qij = / p^x'x^cfx. (59) 

J P>Pcv,t 

Here, the integration is carried out only for p > pcut where pcut is a selected eutoff density. In this paper, we chose 
as Pent = Pmax/100 to focus On the high-density region. For comparison, we also chose the cutoff density as zero (i.e., 
the distortion parameter is defined in terms of 7^ ). In this case, the distortion parameter is denoted as 779. In the 
following, we primarily adopt r], and if the value of r] grows exponentially, we judge that the model is dynamically 
unstable. 

In Fig. 9, we show the evolution of pmax, etc, and 77 for models M5clC, M5c2A, M5c2C, M7b3A, M7b3C, M7c2A, 
M7c2C, M7c3A, and M7c3C. For models M5c2A, M7b3A, M7b3C, M7c2A, M7c3A, and M7c4C, (and also rjo) does 
not increase exponentially (cf. Table IV). This implies that these models are dynamically stable against bar-mode 
deformation. For models M5clC, M5c2C, M7c2C, and M7c3C, on the other hand, the values of t] approximately 
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FIG. 10. Snapshots of the density contour curves for p in the equatorial plane for model M7c2C. The sohd contour curves 
are drawn for p/pmax = 1, 0.8, 0.6, 0.4, 0.2, and 10"-'^^ for j = 2, • • • , 8. Vectors indicate the local velocity field {v^ ,v^), and 
the scale is shown in the upper right-hand corner. 



increase in proportional to e*/"^ where r denotes a characteristic growth time. Thus, these models are unstable. The 
growth time r of these models is approximately 0.8-1 msec. This is the same order of magnitude as the dynamical 

— 1/2 

time scale pc ■ Thus, the instabilities found here are indeed the dynamical instabilities. 

In Figs. 10 12, we display the snapshots of the density contour curves and velocity vectors for the selected time 
slices for models M7c2C, M7c3C, and M5c2C. In all the models, the collapse proceeds in approximately axisymmetric 
manner throughout the initial collapse to the first bounce, forming a torus-like structure. For M7c2C, after the first 
bounce, the formed core expands by a large factor, and then, collapses again. In this second collapse, nonaxisymmetric 
instabilities grow significantly: In the torus-like high-density region, two density peaks are formed (third panel of Fig. 
10). Then, the separation of the density peaks increases, and a bar-like structure is formed (fourth panel), developing 
spiral arms in the outer region. Subsequently, the separation decreases, and they eventually merge and form a single 
peak (fifth and sixth panels). In the outer region, spiral arms are developed, which play a role for transferring the 
angular momentum of the formed core to the outer region. Because of this angular momentum transfer as wc;ll as the 
dynamical friction force to the bar from the surrounding matter, the nonaxisymmetric structure of the central core is 
quickly erased and the protoneutron star eventually relaxes to a slightly nonaxisymmetric quasistationary state. 

For model M7c3C, the nonaxisymmetric instabilities grow in a similar manner to that of M7c2C. In this case, 
however, the maximum value of rj, which denotes the achieved maximum degree of nonaxisymmetric deformation, is 
slightly smaller. This seems to reflect the fact that the angular momentum is not as large as that of M7c2C. For model 
M5c2C, the evolution is very similar to that of M7c3C. However, the growth rate of i] for M5c2C is slightly smaller 
than for M7c2C. The reason is that the mass and the compactness of the outcome formed after the collapse are smaller, 
and hence, the growth time of the nonaxisymmetric dynamical instabilities, which is approximately proportional to 
the dynamical time scale, becomes longer. 

A noteworthy feature for the unstable models is that in the late phase in which the bar-mode perturbation damps, 
the m = I mode grows gradually and becomes a dominate mode eventually. With the growth of this mode, small 
one-armed spiral arm is formed (see, e.g., the last panels of Figs. 10-12). The excitation of this mode is probably due 
to the fact that the formed star is highly differentially rotating [55,56]. However, the effect of its growth is not very 
outstanding since the amplitude of the perturbation is not very large and fairly quickly damps due to the angular 
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momentum transfer to the outer region. Therefore, we conchide that the onset of to = 1 mode instabihties is not as 
important as that of the bar-mode for the evolution of the system, although the density configuration of the formed 
protoneutron star becomes asymmetric due to it. 

The formation of the bar and subsequent outward transfer of the angular momentum change the density profile 
of the protoneutron stars. In Fig. 13, we show the evolution of a.c, Pmax, and t] for models M7c2C and M7c3C 
in the three-dimensional simulations as well as in the axisymmctric ones. To illustrate that the convergence is 
approximately achieved, the three-dimensional results with N = 156, 188, and 220 are shown together. In the early 
stage of the evolution {t < 88 msec) in which the amplitude of the bar-mode perturbation is small, the results of the 
three-dimensional and axisymmetric simulations are in good agreement: i.e., Oc and pmax arc simply in a damped 
oscillation. Slight disagreement between the results of the three-dimensional and axisymmetric simulations is likely 
due to the fact that we discard the matter located in the outer region in the three-dimensional simulations. In a 
stage in which the system is approximately axisymmetric, shock dissipation which damps the oscillation is only the 
mechanism for modifying the density profile. On the other hand, in the late stage with r; > 0.1 (see Figs. 13 (c) and 
(d)), pmax (etc) gradually increases (decreases) with time in the three-dimensional simulations. This reflects the effect 
of the angular momentum transfer by which the centrifugal force in the inner region is weaken and the formed object 
becomes more compact than the outcome in the axisymmetric simulations. In particular, the effect is remarkable for 
model M7c2C. In this case, an oscillating (type 0-A) star (not a protoneutron star in the definition of this paper) is 
formed in the early phase of the axisymmetric simulation, while in the three-dimensional simulation, the protoneutron 
star is promptly formed because of the quick contraction due to the outward transfer of the angular momentum. For 
more massive cases with M ^ 3Mq, protoneutron stars which are supported by strong differential rotation may be 
formed first [59], but the angular momentum transfer may trigger black hole formation. This effect may also play 
an important role in stellar collapse of very massive (M > 25OM0) stars (population III stars) which is triggered 
by the electron-positron pair creation instability [60]. Very massive stars arc likely to be rapidly rotating [61], and 
the collapse may not result directly in a black hole but in very massive self-gravitating disks [62]. The disks will 
be dynamically unstable against nonaxisymmetric deformation and the resulting angular momentum transfer by a 
nonaxisymmetric structure may induce black hole formation. 

Figure 13 also shows that with increasing the value of N, the numerical results achieve a convergence. The results 
oi N = 188 and 220 agree well (except for those in the very late time for which the numerical error is accumulated 
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too much), implying that a convergent result is obtained with A'' ~ 200. By the way, the period of the quasiradial 
oscillation becomes spuriously longer due to the larger numerical dissipation with the smaller value of N. Such 
spurious effect may lead to underestimation of the growth rate and the achieved maximum value of 77. The convergent 
test carried out here gives us a caution that we have to guarantcx; a sufficient grid resolution in this problem. 

In the lower panels of Fig. 13(c) and (d), we display the evolution of 770 to compare with that of 77. For model 
M7c2C, rjQ increases to be much larger than the initial value. However, the value of rjo is smaller than that of rj even 
when the bar-mode grows to a nonlinear regime. This implies that the bar structure is formed mainly in the central 
region. This feature is more outstanding for model M7c3C in which the increase of rjo from the initial value is not 
seen. Thus, we conclude that the bar-mode perturbation is amplified only in the central region. This is reasonable 
since in the models with ^4 = 0.1, the outcomes are rapidly rotating only in the central region. 



B. Criterion for the onset of nonaxisymmetric dynamical instabilities 

Models M7c2C and M7c3C are dynamically unstable against bar-mode and m = 1 mode deformation, while model 
M7c4C is stable for both modes. This implies that for the onset of dynamical nonaxisymmetric instabilities, high 
values of /3 are necessary for given values of Fi and M. 

All the models with Fi = 1.3 that we picked up here are stable. On the other hand, the models with Fi = 1.28 are 
much more prone to be unstable. This suggests that only for a sufBcicntly small value of Fi < 1.28, the collapsing star 
can be unstable. As mentioned in Sec. IV, for the smaller value of Fi (for the larger depletion factor of the internal 
energy and the pressure at the onset of collapse), the outcomes are more torus- like than those for other values of Fi, 
and also, the degree of differential rotation is larger. These facts are likely reasons that models with Fi = 1.28 are 
more subject to the dynamical instabilities. 

Even with M w 1.5Mq (models M5clC and M5c2C), nonaxisymmetric bar-mode instabilities set in, although 
the maximum values of /? and compactness for the outcomes are smaller than for models M7c2C and M7c3C. This 
indicates that the mass and compactness achieved in the collapse are not very important parameters for triggering 
the bar-mode instabilities as far as M is larger than ^ 1.5Mq. However, it should be noted that general rclativistic 
effects certainly help making a compact outcome. Thus, if the mass is much smaller than ^ I.HMq, nonaxisymmetric 
instabilities may not set in. 

Although a high value of P is necessary, the onset of the nonaxisymmetric dynamical instabilities is not simply 
determined by the value of (3; i.e., although a large value of /3 > 0.27 is preferable for the onset, it is neither the 
necessary nor the sufficient conditions. The first evidence for this statement is that the values of rj for the unstable 
models do not increase at the first bounce at which the value of /? becomes maximum with /3max ^ 0.27. The growth 
of the perturbation is significantly induced in the subsequent bounce stages. Also, model M7c2A is dynamically stable 
although /3niax ~ 0.33. These show that even if (3 exceeds ^ 0.27, the nonaxisymmetric perturbations do not grow. 
Probably, the duration of the phase for which /3 > 0.27 would have to be much longer than the dynamical time scale 
for the onset of the dynamical instabilities. 

Second, the value of /? during the growth of the bar-mode perturbation is smaller than 0.27 for any unstable model. 
In Fig. 14, we show the time evolution of (3 for models M7c2C, M7c3A, M7c3C, and M5c2C. It shows that during 
the growth of the perturbation, f3 for models M7c3C and M5c2C is at most 0.25 and in average 0.2, which is 
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FIG. 13. Evolution of and pmax (a) for model M7c2C and (b) for model M7c3C in three-dimensional simulations with 
A'' = 156 (dashed curves), 188 (solid curves), and 220 (long-dashed curves) as well as in axisymmetric simulation (dotted 
curves). Evolution of rj and rjo (c) for model M7c2C and (d) for model M7c3C in three-dimensional simulations with varying 
grid resolution. 



much smaller than the widely-believed critical value ~ 0.27. There are at least three possible reasons that /3 may be 
smaller than 0.27 for the onset of the nonaxisymmetric dynamical instabilities. The first one is that the effective value 
of /3 in the high-density region may be larger than the global value, and may be large enough for the onset of the 
nonaxisymmetric dynamical instabilities. This is likely to be the case in particular for the highly differentially rotating 
collapse since the effective value of (3 in the central high-density region where the nonaxisymmetric perturbation 
grows dominantly is larger than the whole value for such cases. The second possibility is that the onset of the 
nonaxisymmetric dynamical instabilities is due to the high degree of differential rotation as indicated in a Newtonian 
simulation [28]. In this case, the high value oi (3 > 0.27 is not necessary. The third possibility is that general 
relativistic effects reduce the critical value of f3 below 0.27. Indeed, in [15], we showed that the critical value of /3 can 
be decreased by ~ 10% due to the general relativistic effects for compact stars with the compactness ~ 0.1-0.2. All 
these possibilities show that the critical value of P for the onset of the nonaxisymmetric dynamical instabilities may 
be smaller than 0.27 depending sensitively on several parameters, and thus, it cannot be uniquely determined. 

Figure 14 also shows that the values of (3 for models M7c3A and M7c3C are not very different during the oscillation 
phase although they are stable and unstable against bar-mode deformation, respectively. This also illustrates that the 
value of f3 does not uniquely determine the dynamical stability. As shown in Fig. 7, on the other hand, the profiles of 
the density and the rotational angular velocity in the central region are different between two models. Thus, in this 
case, the degree of differential rotation and the steepness of the density profile play an important role for determining 
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FIG. 14. Evolution of /3 for models M7c2C (solid curve), M7c3C (long-dashed curve), M7c3A (dashed curve), and M5c2C 
(dotted-dashed curves) in three-dimensional simulations. The values of t2d are listed in the caption of Fig. 9. The dotted line 
denotes f3 = 0.27. 



the stability. 

The case with Fi = 1.28, in which the depletion of the pressure and the internal energy in an early stage of collapse 
with p <C Pnuc is largest among the three cases, is more subject to the nonaxisymmetric dynamical instabilities. This 
indicates that a large depletion of the internal energy and the pressure in the early stage is an essential element 
for the onset of the nonaxisymmetric dynamical instabilities. The reason is that for the larger depletion factor, the 
collapse in the central region proceeds significantly to make a compact core, and hence, to increase the spin of the 
central region as illustrated in the axisymmetric simulations (cf. Fig. 7). In a realistic phenomena, the depletion of 
the pressure and the internal energy in the early stage is determined by the partial photo-dissociation of the iron to 
lighter elements, by the electron capture, and by the neutronization [48,49]. Since the depletion factor is a crucial 
parameter, an appropriate modeling for such microphysical processes will be necessary for a more detailed study on 
the nonaxisymmetric dynamical instabilities in the future. 

No evidence for fragmentation of protoneutron stars is found in the first ^ 10-20 msec after the bounce in the 
present numerical simulations. Previous studies in the field of protostar formation from collapsing gas clouds (e.g., 
[63]) show that the fragmentation occurs when the thermal energy at an initial stage of the collapse is much smaller 
than the gravitational potential energy: In the case of the small thermal energy, a torus-like or a disk-like structure is 
formed as a result of the collapse and subsequently the fragmentation takes place. This indicates that if the value of 
Fi is much smaller than 1.28 (i.e., if the fraction of the depletion of the internal energy and the pressure in an early 
stage of collapse is much larger than 16%), the fragmentation may occur during the stellar core collapse. However, 
such an extremely small value of Fi (an extremely large value of the depletion factor) is unlikely to be achieved in 
the stellar core collapse [48,49], and therefore, we infer that the fragmentation of protoneutron stars would not occur 
in nature, at least in a few 10 msec after the stellar collapse. 

To summarize, the nonaxisymmetric dynamical instabilities set in only for the case that the following conditions 
are satisfied: (i) the progenitor of the stellar core collapse is rapidly rotating with the initial value of /3 > 0.01, (ii) 
the degree of differential rotation for the velocity profile of the initial condition is very high with A < 0.1, (iii) the 
depletion factor of the pressure and the internal energy in an early stage of collapse in which p <C Pnuc should be 
large enough to induce a rapid collapse in the central region of the stellar core and for an efficient spin-up. With the 
increase of stellar core mass, the maximum value of f3 achieved during the collapse is increased, but this does not 
significantly change the stability property as far as M is larger than ~ 1.5Mq. It is also found that the value of /3 
does not uniquely determine the property of the dynamical stabilities. 

C. Gravitational waveforms from nonaxisymmetrically deformed stars 

In Fig. 15, we show gravitational waveforms and total emitted energy and angular momentum as a function of 
retarded time for models M7c2C, M7c3C, and M5c2C. For these models, nonaxisymmetric dynamical instabilities set 
in after the bounce resulting in formation of a bar and spiral arms and in excitation of gravitational waves with m = 2 
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FIG. 15. Gravitational waveforms in the gauge-invariant wave extraction method (i?+,x) and in the quadrupole formula 
(A+,x) (a) for M7c2C, (b) for M7c3C, and (c) for M5c2C. (d) Total emitted energy and angular momentum as a function of 
retarded time for models M7c2C (solid curves), M7c3C (dashed curves), and M5c2C (long-dashed curves). 



modes. Gravitational waveforms are computed both by the gauge-invariant wave extraction and by the quadrupole 
formula. 

Figures show that with the amphfication of r], the amplitude of gravitational waves are increased. However, once it 
reaches the maximum, the amplitude damps quickly as in the evolution of r]. This is due to the effect that the bar-mode 
perturbation plays a role for transferring the angular momentum from the inner region to the outer one. Eventually, 
the bar-mode perturbation damps, resulting in a quick damping of gravitational wave amplitude. The damping is 
in particular outstanding for model M7c2C. This is due to the fact in this model, the amplitude of the bar-mode is 
largest, and hence, the angular momentum transfer is most effective. Due to this fact, the total emitted energy for 
models M7c2C and M7c3C becomes approximately identical although the maximum amplitude of gravitational waves 
for M7c2C is about twice larger. 

In isolated rotating stars, once the bar-mode instabilities set in and saturate, the amplitude of their perturbation 
remains approximately constant, resulting in emission of quasipcriodic gravitational waves in a dissipation time scale 
of gravitational radiation which is much longer than the dynamical time scale (e.g., [28]). However, in the rotating 
core collapse, the amplitude of gravitational waves is damped by the angular momentum transfer from the bar to the 
surrounding matter, for which the time scale is nearly equal to the dynamical time scale and much shorter than the 
emission time scale of gravitational radiation. 

The maximum amplitude of gravitational waves for model M7c2C is by a factor of ~ 2 larger than that for M7c3C, 
reflecting that the degree of nonaxisymmetric deformation is larger. The amplitude for model M5c2C is by a factor 
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of ^ 2.5 smaller than that of M7c3C, although the initial value of (3 is approximately identical and the waveforms 
are very similar for these two models. According to the quadrupole formula, the amplitude of gravitational waves is 
approximately proportional to if the radius of the formed protoneutron star is identical. Thus, the dependence 
on mass is reflected in the amplitude. 

In the evolution of models M7c2C, M7c3C, and M5c2C, the m = 1 mode perturbation grows in the late phase of 
the evolution. However, this does not affect the amplitude of gravitational waves significantly, since the amplitude of 
the perturbation is not very large and the m = 1 mode does not contribute to the lowest-order (mass quadrupole) 
waveforms in the three-space of the reflection symmetry with respect to the equatorial plane. 

For models M7c2C and M7c3C, the maximum values of i?+,x are ~ 0.25 km and 0.15 km, respectively. For M5c2C, 
it is even smaller ~ 0.05 km. The amplitude of gravitational waves, h, observed at a distance of r along the optimistic 
direction {6 = 0) is written as 



21/^ ^+,x W lOMpc \ 
Vo.aikmyV r J' 



(60) 



This implies that the observed amplitude at a distance of 10 Mpc is at most /i < 8 x 10 for initial core mass 
M ~ 2.5M0 and /i ~ 1.5 x IQ-^^ for M ~ 1.5Mq. 
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FIG. 16. Ao in the axisymmetric simulations (a) for models M7al (solid curves) and M5al (dashed curves), and (b) for 
models M7c3 (solid curve) and M5c2 (dashed curve) with Fi = 1.3 and F2 = 2.5. 



To compare the amplitude of gravitational waves from the bar-mode deformation with that from axisymmetric 
collapse, we show Aq for models M5al, M7al, M5c2, and M7c3 with Fi = 1.3 and r2 = 2.5 in the axisymmetric 
simulations in Fig. 16. As mentioned in [12], it is difficult to extract gravitational waves of small amplitude from the 
metric in the axisymmetric simulations, and hence, only the waveforms by the quadrupole formula are presented here. 
Although it provides only an approximate waveform, the wave phase can be accurately computed and the error of the 
amplitude will be at most ^ 10%. Figure 16 indicates that for the initial mass M ^ 2.5Af0, the maximum amplitude 
is at most 0.01 km for the rigidly rotating case and 0.02-0.03 km for differentially rotating cases. The values are by 
a factor of ^ 2 smaller for M ~ 1.5Mq. Thus, the amplitude of gravitational waves of the I = m = 2 modes from 
the nonaxisymmetric dynamical instabilities is ^ 10 times as large as that in the axisymmetric case. On the other 
hand, those amplitudes are not as large as the maximum amplitude of gravitational waves from of coalescing binary 
neutron stars in close circular orbits [36]. Thus, the nonaxisymmetric deformation in the stellar core collapse is not 
as a stronger emitter as coalescing binary neutron stars. 

For model M7c2C (M7c3C), the total emitted energy and angular momentum are about 9 x 10^° erg (9 x 10^° erg) 
and 3 x 10''^ g cm^/sec (2 x 10^^ g cm^/scc), respectively. These values are about 0.03% (0.03%) of the total mass 
energy (M*c^) and 0.7% (0.6%) of the total angular momentum, respectively, and are much larger than those in the 
axisymmetric collapse (e.g., [10]). However, they are not as large as those in merger of binary compact objects in 
which ^ 1% of the total mass energy and 10% of the angular momentum are dissipated by gravitational waves 
in the final phase of the merger [36]. Thus, in the stellar collapse, the radiation reaction by gravitational waves are 
not likely to play an important role for the dynamics of bounce and oscillation of the protoneutron star. For model 
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M5c2C, these values are much smaller because of its small mass and small compactness achieved. Hence, the effect 
of gravitational wave emission is less important. 
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FIG. 17. Energy power spectra oi I — m — 2 modes (a) for models M7c2C (solid curve) and M7c3C (dashed curve), and (b) 
for M5c2C. 
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FIG. 18. Gravitational waveforms and energy luminosity (a) for model M7c2C and (b) for model M7c3C with A'^ 
and 220. The long-dashed, solid, and dashed curves denote the results with A'^ = 220, 188, and 156. 
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Comparison of gravitational waveforms computed by the gauge-invariant wave extraction method and by the 
quadrupole formula shows that the wave phase in the two results agree approximately (besides a systematic phase 
shift). However, the amplitude disagrees by a factor of ^ 2. As pointed out in [12], in the quadrupole formula, the 
amplitude is underestimated by a factor of Al/R ^ 0.1-0.2 where R here denotes the characteristic radius of the 
outcome after the collapse. On the other hand, the amplitude in the gauge-invariant wave extraction method is likely 
to be overestimated because the waveforms are extracted in a local wave zone [44]: In this paper, L ^ A/2 < A where 
A is the wavelength of gravitational waves ^ 300 km, and thus, the amplitude would be overestimated by a factor 
of 10-20% [44]. Hence, the true amplitude would be between two results. However, besides the disagreement in the 
amplitude, two methods provide qualitatively the identical results. This reconfirms that the quadrupole formula is a 
reasonable method for approximately computing gravitational waveforms even in fully general relativistic simulations, 
in the absence of black holes. 
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In Fig. 17, wc show the energy power spectrum of m = 2 modes (a) for models M7c2C (solid curve) and M7c3C 
(dashed curve) and (b) for model M5c2C. The Fourier spectrum is computed from gravitational waveforms in terms of 
the gauge-invariant wave extraction. The spectrum for model M7c2C is broader in a low frequency region with / < 1 
kHz than those for other models. This reflects the long oscillation period of this model. The peak frequency is about 
0.8-1.3 kHz in these models. These frequencies are determined by the quadrupole f mode frequency of the deformed 
star formed after the bo unce. Namely, the higher peak frequency implies that the outcome is more compact in 
proportional to \JM/R^. According to a perturbative study for the quadrupole f mode [64] , the frequency of neutron 
stars becomes ~ 2.5-4 kHz. The frequency of the oscillation of unstable protoneutron stars is much lower than that 
of neutron stars. The reason is that the radius of the protoneutron star is larger. Nevertheless, the peak frequency 
is higher than the best sensitive frequency (between ~ 100 and several 100 Hz) of kilometer size laser-interferometers 
such as LIGO [34]. As shown in Eq. (60), the amplitude of gravitational waves is not very high if we assume the 
distance to the source ^ 10 Mpc. Thus, gravitational waves from nonaxisymmetrically deformed protoneutron stars 
may be promising sources for such gravitational wave detectors only when the stellar collapse happens for r <C 10 
Mpc. On the other hand, the frequency may be in a good range for resonant-mass detectors and/or specially designed 
advanced interferometers such as the advanced LIGO [34]. 

To summarize this section, we have found that the amplitude of gravitational waves from dynamically unstable 
protoneutron stars against nonaxisymmetric deformation is ^ 10 times as large as that from the axisymmetric collapse. 
However, even in the case that the degree of the nonaxisymmetric deformation is as large as in model M7c2C, the 
maximum amplitude is < 20-30% of that in merger of binary neutron stars (e.g., [36]). Since the peak frequency 
of gravitational waves is fairly high ^ 1 kHz, gravitational waves from nonaxisymmetric dynamical deformation of 
protoneutron stars may become promising sources for the laser-interferometric gravitational wave detectors only in 
the case that the event rate for the nonaxisymmetric deformation in the stellar core collapse is large. 

Before closing this section, we demonstrate that the convergence with improvement of the grid resolution is achieved 
fairly well for gravitational waveforms. In Fig. 18, we show the numerical results for models M7c2C and M7c3C with 
N = 156, 188, and 220. For the lower grid resolution, the period of the quasiradial oscillation becomes longer. As a 
result, the growth rate of rj becomes smaller. This causes an error in phase of gravitational waves. Also, the lower 
grid resolution results in underestimating the maximum value of r}. As a result, the amplitude of gravitational waves 
is underestimated. However, with N > 200, the numerical results appear to converge well. Thus, we conclude that 
with our choice of the grid resolution, a good convergent result is obtained. 

VI. SUMMARY AND DISCUSSION 

We have presented the first numerical results of three-dimensional hydrodynamic simulations for stellar core collapse 
in full general relativity focusing mainly on the criterion for the onset of the bar-mode dynamical instabilities. The 
nonaxisymmetric dynamical instabilities have been widely studied for isolated rotating stars in equilibrium to this time 
not only in Newtonian gravity but also in general relativity. However, for nonaxisymmetric dynamical instabilities in 
rotating stellar core collapse, very little study has been done even in Newtonian gravity [7] . Taking into account such 
status, we performed the simulations for a wide variety of equations of state, stellar masses, and velocity profiles to 
clarify the criterion for the onset of the nonaxisymmetric dynamical instabilities as well as the outcomes after their 
onset. 

A number of previous works for isolated rotating stars in equilibrium have clarified that the bar-mode dynamical 
instabilities can set in when the value of (i exceeds ~ 0.27 or when the degree of differential rotation is sufficiently 
high. Thus, first, we performed axisymmetric simulations of rotating stellar collapse to clarify the conditions that 
the value of (3 is amplified beyond ^ 0.27 and that the degree of differential rotation for the outcomes of the collapse 
becomes very large. We have found the following conditions are necessary to achieve a state with /3max > 0.27: (A) 
the initial state of the collapse is highly differentially rotating with A < 0.1; (B) the progenitor is moderately rapidly 
rotating with 0.01 ^ /?i„it ^ 0.02, but has to be not very rapidly rotating such as /3i„it ^ 0.02; (C) the progenitor 
star is massive enough to achieve a compact state for which a significant spin-up is achieved. However, at the same 
time, the mass should not be very high to avoid black hole formation. We also found that to achieve a high degree of 
differential rotation after the collapse, the depletion factor of the pressure and the internal energy in an early stage of 
collapse in which p <C Pnuc should be large enough to induce a rapid collapse in the central region of the stellar core 
and resulting efficient spin-up. 

Next, to clarify the condition for the onset of nonaxisymmetric dynamical instabilities, we performed three- 
dimensional simulations. Based on the results of axisymmetric simulations, we picked up models which are likely 
to become unstable during the collapse and bounce. From the three-dimensional simulations, it is found that the 
nonaxisymmetric dynamical instabilities set in only for a restricted parameter range as indicated by axisymmetric 
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simulations. Specifically, the following conditions are required to be satisfied: (i) the progenitor of the stellar core 
collapse is rapidly rotating with 0.01 /3init ^ 0.02, (ii) the degree of differential rotation for the velocity profile of 
the initial condition is very high with A ^ 0.1, and (iii) the depletion factor of the pressure and internal energy in the 
early stage of collapse is large enough to induce a rapid collapse in the central region of the stellar core. Although 
stellar cores of larger mass have more advantage to form a compact protoneutron star, the stability property depends 
weakly on the mass as far as M > I.SMq. 

Gravitational waves are computed in the case that the bar-mode dynamical instabilities set in. For the case that 
the bar-mode perturbation grows, the amplitude of gravitational waves increases exponentially, and as a result, burst- 
type waves arc emitted. However, since the bar-mode of the core subsequently damps due to the outward angular 
momentum transfer in a short time scale ~ 2-3 msec, the amplitude of gravitational waves decreases quickly. Thus, 
quasiperiodic gravitational waves, in which the amplitude can be accumulated effectively, are not emitted efficiently 
after the damping of the nonaxisymmetric perturbation. The maximum amplitude of gravitational waves at a distance 
of 10 Mpc is ~ 4-8x10"^^ with the frequency ^ 1 kHz for very massive core collapse with initial core mass M ~ 2.5Mq. 
The maximum amplitude is approximately proportional to for a given value of Anit- For M ~ 1.5Mq, thus, the 
maximum amplitude is ~ 1-2 x 10~^^ at a distance of 10 Mpc. This amplitude is about 10 times as large as that 
in the axisymmetric collapse, but ~ 20-30% of the maximum amplitude in the merger of binary neutron stars (e.g., 
[36]). Thus, the feature of gravitational waves is summarized as follows: (i) burst-type (not quasiperiodic) waves are 
emitted, (ii) the frequency is relatively high with ~ 1 kHz, and (iii) the amplitude is about 10 times as large as that 
from axisymmetric collapse, but not as large as that for merge; of binary compact objects. These facts imply that 
only when nonaxisymmetric dynamical instabilities set in for a large fraction of the stellar core collapse, gravitational 
waves induced by nonaxisymmetric dynamical instabilities of protoneutron stars may become promising sources for 
kilometer size laser-interferometers. 

Besides the dynamical instabilities, there is another route for nonaxisymmetric deformation: secular instabilities. 
As found from Figs. 1 and 6, the value of f3 in the protoneutron stars formed after the collapse is often larger than 
the critical value for the onset of secular instabilities ^ 0.14. According to previous works [65-67], isolated rotating 
stars of /3 > 0.14 can form an ellipsoidal structure due to gravitational radiation, which may become a strong emitter 
of quasiperiodic gravitational waves with the frequency between 10 and several 100 Hz. However, these studies were 
performed for isolated stars. In the case of stellar core collapse, the formed protoneutron stars will be surrounded 
by massive outer envelopes, and thus, the bar-mode perturbation excited by the secular instabilities may be damped 
quickly due to the angidar momentum transfer from the ellipsoidal protoneutron star to the outer envelope as in 
the case of the dynamical instabilities. A simulation with massive envelope will be necessary to clarify whether the 
secular instabilities can grow or not. On the other hand, in contrast to the dynamical instabilities, the growth time 
scale of the secular instabilities is fairly long > 100 msec. In such a long time scale, the surrounding matter may be 
ejected outward or accreted onto the central neutron star in reality, and hence, the secular instabilities may grow as 
in the isolated stars. However, in such a long time scale, viscous or magnetic dissipation may also play an important 
role [59] for preventing the growth of the nonaxisymmetric perturbation. At present, it is totally unclear whether the 
secular instabilities set in or not. 

As reported in this paper, the criterion for the onset of nonaxisymmetric dynamical instabilities may depend sen- 
sitively on the equations of state for subnuclear density, since with the smaller pressure for p < pnuc, the collapse is 
accelerated more for an efficient spin-up of the central region. In the present work, we adopted a parametric equation 
of state for simplicity. To clarify the criterion for the onset of nonaxisymmetric dynamical instabilities more strictly, 
however, sophisticated equations of state should be adopted. In realistic equations of state, the increase rate of the 
pressure as a function of the density (i.e., an adiabatic index) significantly decreases at density ~ lO^^g/cm^ (e.g., 
[68,69]). This suggests that in a realistic equation of state, collapse of the central region is likely to be accelerated sig- 
nificantly before reaching the nuclear density and, hence, the collapsed core may be more subject to nonaxisymmetric 
dynamical instabilities. This fact suggests that a simulation with more realistic equations of state is an interesting 
subject for the future. 

Finally, we note the following issue. This paper focuses only on nonaxisymmetric dynamical instabilities of pro- 
toneutron stars in the first 10-20 msec after the bounce. The formed protoneutron stars subsequently emit neutrinos 
and dissipate the thermal energy [49,70]. As a result, they contract gradually in a time scale of ~ 10 sec. Because 
the angular momentum is conserved approximately, the spin of the protoneutron stars will be increased with the 
contraction and (3 may be increased beyond ^ 0.27. Thus, even in the case that they are stable in the first 10-20 
msec after the bounce, they may eventually become unstable after the neutrino cooling. This issue is not investigated 
in this paper and remains for the future. 
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